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ABSTRACT 

A new algorithm for identifying and characterizing vortices in complex flows is presented. 
The scheme uses both the vorticity and pressure fields. A skeleton line along the center of a 
vortex is produced by a two-step predictor-corrector scheme. The technique uses the vector 
field to move in the direction of the skeleton line and the scalar field to correct the location in 
the plane perpendicular to the skeleton line. A general vortex cross section can be concisely 
defined with five parameters at each point along the skeleton line. The details of the method 
and examples of its use are discussed. 
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1 Introduction 


Vortices are considered the most important structures that control the dynamics of flow 
fields. Large-scale vortices observed in atmospheric and oceanographic flows are responsible 
for ozone holes, hurricanes, tornadoes, and maelstroms. Vortices that are shed from aircraft 
canards, wings, and control surfaces affect the handling characteristics of the airplane and 
determine the closest safe following distance of other aircraft. On a smaller scale, vortices are 
the fundamental building blocks of turbulent flow structures and are intimately connected 
with the mixing of chemical species, heat transfer, and drag forces. 

Although the term “vortex” connotes a similar concept in the minds of most fluid dynam- 
icists, a precise definition is difficult to obtain. Robinson [1] suggests the following working 
definition: “A vortex exists when instantaneous streamlines mapped onto a plane normal to 
the vortex core exhibit a roughly circular or spiral pattern, when viewed from a reference 
frame moving with the center of the vortex core.” Unfortunately, this definition requires a 
knowledge of the vortex core before one can determine whether something is a vortex. In 
practice, Robinson [2] and Robinson, Kline, and Spalart [3] use the above rigorous definition 
to confirm that a particular structure is, in fact, a vortex. Regions of low pressure are used 
to identify candidate vortices. Their scheme exploits the fact that the pressure in the core of 
a vortex is lower than that of the surrounding fluid. A radial force is needed to provide the 
centripetal acceleration that keeps a particle rotating about an axis. In the case of vortices, 
this force is in the form of the pressure; the pressure inside a vortex is lower than the pressure 
outside the vortex. Robinson and his colleagues find that elongated low-pressure regions in 
incompressible turbulent flows almost always indicate vortex cores, but these surfaces can 
be difficult to characterize and provide no information in regard to the sense of rotation or 
the connectivity between the structures. 

Moin and Kim [4] and Kim and Moin [5] use vorticity lines (sometimes called vorticity 
field lines or vortex lines) to visualize vortical structures in turbulent channel flow. A vorticity 
line is everywhere parallel to the local vorticity vector and is defined by 


dx c o 

ds |u3| 


(i) 


where Co — V x u is the vorticity vector, u is the velocity vector, x is the location in space, and 
s is the distance along the vorticity line. The definition suggests that vorticity lines should 
go through the cores of vortices, which is usually the case. However, the resulting vortex 
lines are extremely sensitive to the choice of initial location for the integration. As Moin and 
Kim [4] point out, “If we choose x 0 arbitrarily, the resulting vortex line is likely to wander 
over the whole flow field like a badly tangled fishing line, and it would be very difficult to 



identify the organized structures (if any) through which the line may have passed.” Kim and 
Moin [5] illustrate the potential tangle in their Fig. 2. To avoid such a confusing jumble, they 
carefully select the initial points (inside suspected vortices) and focus on specific parts of 
the resultant lines to illustrate some important features of turbulent channel flow. However, 
Robinson [1] shows that even experienced researchers can be surprisingly misled by ordinary 
vorticity lines. 

Villasenor and Vincent [6] present an algorithm for the recognition and visualization of 
elongated vortices in three-dimensional time-dependent flow fields. They start from a seed 
point and compute the average length of all vorticity vectors contained in a small-radius 
cylinder. The axis of the cylinder extends from the seed point to the surface of a sphere 
centered about the seed point. They repeat this step for a large number of cylinders that 
emanate from the seed point. The cylinder with the maximum average vorticity vector length 
becomes a segment of the vortex tube. This process is repeated until the vorticity decreases 
enough to end the vortex tube. Because the scheme depends only on the magnitude of 
the vorticity vectors rather than their direction, the algorithm must limit the angle between 
adjacent line segments so that the new line segment is different from the previous line segment 
traced in the opposite direction. All vorticity vectors inside the vortex tube are drawn. The 
process can be repeated for multiple time steps so that the evolution of the vortex tubes is 
visualized as a sequence of all stored vorticity vectors. Villasenor and Vincent [6] claim that 
the use of only the magnitudes and not the directions of the vorticity vectors is advantageous 
because it allows the algorithm to find structures that do not have axes aligned with the 
principal axis of the structure. As a consequence, the algorithm can inadvertently capture 
structures that are not vortices. In addition, by visualizing all vectors that are inside the 
cylinder, structures are included that are not part of the vortex originally intended to be 
viewed. An example can be seen in Fig. 4(c) of Villasenor and Vincent [6], in which another 
structure is responsible for some of the vorticity vectors that are almost orthogonal to the 
vortex axis. 

Zabusky, Boratav, Pelz, Gao, Silver, and Cooper [7] fit ellipsoids to surfaces of constant 
vorticity |u;| and to constant vortex stretching |u ; • Vt?|/|u;| in an effort to understand the 
dynamics of a vortex reconnection process. A bundle of vector field lines of the vorticity 
u; and the vortex stretching tu ■ Vu emanate from the major and minor axes and from the 
center of the minor ellipse of the respective fitted ellipsoids. Both the ellipsoids and the 
vector field lines provide useful information about the vortices for their flow field, which 
contains neither solid boundaries nor a mean straining field. In flows with these additional 
complications, the regions with large vorticity magnitudes do not necessarily correspond to 
vortices; hence, ellipsoids of constant vorticity coupled with vorticity lines are unlikely to 
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provide useful characterizations. 

Chong, Perry, and Cantwell [8] address the question of when a region of vorticity is 
a vortex. They suggest that a vortex core is a region where the velocity-gradient tensor 
has complex eigenvalues. In such a region, the rotation tensor dominates over the rate-of- 
strain tensor. Soria and Cantwell [9] use this approach to study vortical structures in free- 
shear flows. At those poi 11 is at which the vorticity magnitude exceeds a specified fraction 
of the maximum vorticity, the eigenvalues of the velocity-gradient matrix are determined. 
Those regions that have complex eigenvalues and satisfy the vorticity-lhreshold criterion are 
rendered as solid surfaces that represent the vortices in the flow. This method correctly 
identifies the large vortical structures in the flow. However, the method clutters the picture 
with many smaller structures and does not provide a way to link the smaller vortical volumes 
with the larger coherent vortices of which they might be a part. 

Yates and Chapman [10] carefully explore two definitions of vortex cores. One definition 
associates a vortex core with a streamline of minimum curvature within a region of spiraling 
streamlines. The other definition regards the vortex core as the line defined by the local 
maxima of normalized helicity. Under certain circumstances, both definitions produce the 
same vortex core. Unfortunately, all of the analyses and conclusions are appropriate only 
for steady flows. 

Bernard, Thomas, and Handler [11] use a semiautomated procedure to identify quasi- 
streamwise vortices in the near-wall region of turbulent-channel-flow direct-numerical-simulation 
data. Their method finds local centers of rotation in user-specified regions in planes perpen- 
dicular to the streamwise direction. The local centers of rotation are linked to corresponding 
centers in adjacent streamwise-normal planes to form a vortex-core line. The results of their 
work indicate that experienced users can correctly find essentially all critical vortices re- 
sponsible for the maintenance of the Reynolds stress. However, their method fails to capture 
vortices with axes that are not nearly aligned with the streamwise direction. This shortcom- 
ing can be a major difficulty for flows such as free-shear layers and transitional boundary 
layers, in which both streamwise- and spanwise-oriented vortices are important. In addition, 
the procedure depends heavily on user intervention to indicate regions where the program 
should search for centers of rotation. This process is tedious, and the detailed results depend 
upon nonquantifiable whims of the user. 

Here, we present a new approach for identifying and characterizing vortices in complex 
flow fields. Rather than a dependence upon a single quantity to define the vortex, the 
new approach is a predictor-corrector method that converges to the vortex skeletons (the 
lines that pass through the vortex cores). In section 2, we describe how the new method 
determines the vortex skeletons. In section 3, we discuss some novel vortex cross-section 
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schemes that put meat on the skeleton lines, and we summarize our findings in section 4. 


2 The vorticity-predictor pressure-corrector method 

Our new predictor-corrector method produces an ordered set of points that approximates 
the skeleton of a vortex. Associated with each point are quantities that describe* the local 
characteristics of the vortex. These quantities may include the vorticity, the pressure, cross- 
sectional information, the circulation, and other quantities of interest. 

The new method produces lines that are similar to vorticity lines. However, because a 
distinction exists between the vorticity, which is a mathematical function of the instantaneous 
velocity field, and a vortex, which is a physical structure witfi coherence over a region of 
space, the straightforward integration of the vorticity lines must be modified to locate the 
skeleton lines associated with individual vortices. At the very least, because vorticity lines 
can begin and end only at domain boundaries and vortices have no such restriction, we must 
develop criteria for initiating and terminating the integration. In addition, the method must 
be self-correcting (i.e., line trajectories that diverge from the vortex core should be subject 
to a restoring (correcting) force). 

In this section, we will first discuss the procedure used to find an initial point on the 
vortex-core skeleton. Then we will reveal the details involved in growing the vortex-core 
skeleton from the seed point. Finally, the termination of the vortex skeleton will be addressed. 

2,1 Finding a seed point 

Any flow-field information can be used to find starting locations for the predictor-corrector 
process. The pressure and the vorticity magnitude are convenient scalars for choosing seed 
points. In our work, the flow field is scanned in planes perpendicular to the streamwise 
direction. The scanning direction aiTects the order in which vortices are located, but not the 
overall features of the vortices. In each plane, the values of the pressure and the vorticity 
magnitude are checked. Threshold values can be chosen a priori, or they can be a predeter- 
mined fraction of the extrema. If the pressure is less than the pressure threshold and the 
vorticity magnitude is greater than its threshold, then the point location is further refined. 
In our implementation, this positional refinement shifts the seed point to the location of the 
local pressure minimum in the plane perpendicular to the vorticity vector. (The rationale 
for this refinement will be discussed later.) The sampling can be done at every grid point. 
We experimented with sampling at every second, third, and fourth points. Although a test 
of each grid point to determine if the specified pressure and vorticity criteria are met is not 
expensive, we had hoped that sampling for seed points on a coarser mesh would eliminate 
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the multitude of redundant vortices that we found. Global coarsening of the grid does not 
uniformly eliminate the problem; some of the redundancies are avoided, but others remain. 

A better method for eliminating the redundant vortices involves the selective elimination 
of grid points from the supply of possible seeds. Elements of a three-dimensional integer 
array correspond to volumes about each grid point. All of the elements are originally set 
to zero. For the first vortex, any point in the flow field may be evaluated as a potential 
seed point. As points are added to the vortex skeleton, the elements of the integer array 
associated with the points are (lagged. The flagged points will not be available as seeds 
for subsequent vortices. Although many coincident vortices are eliminated in this way, a 
further reduction can be achieved by also flagging elements of the array that correspond to 
points that are within the vortex cross section. In Fig. 1, multiple skeleton lines started from 
different seed points all pass through the same vortex tube. The redundancies are eliminated 
when points inside a tube are excluded from the pool of future seed points. 

2.2 Growing the skeleton 

Once a seed point has been selected, the skeleton of the vortex core can be grown from 
this point. The scheme that we have developed is a two-stage predictor-corrector method. 
With this technique, the next position of the vortex skeleton is predicted by integrating 
along the vorticity vector. This candidate location is corrected by adjusting the position 
to the pressure minimum in the plane that is perpendicular to the vorticity vector. To 
ensure that the minimum is actually part of the vortex under consideration, only a limited 
angle between the vorticity vectors at the predicted and the corrected point is allowed. 1 
The continuous modification of the vorticity-line location lessens the sensitivity to both the 
initial conditions and the integration details. These sensitivities are common problems in 
vorticity-line calculations. Just as importantly, the corrector stage decreases the influence 
of background mean shear and small-scale vorticity fluctuations on the identification of the 
vortex cores. The identified vortex skeleton corresponds more closely to what one might 
intuitively choose as the vortex skeleton than to the standard vorticity-line result. 

The modified vortex-line scheme is illustrated in the schematic diagrams of Fig. 2 and is 
also summarized in pseudocode below. 


l This differs from the angle limitation imposed by Villasenor and Vincent [6]. They restrict the change 
in direction of the skeleton line along the center of the vortex tube from one point on the skeleton line to 
the next. We restrict the change in the direction of the vorticity vector from the predicted location to the 
corrected location. 


For each remaining seed point 

If the seed point p 0 is not in any previous vortex 

While the vortex skeleton continues 

1. Determine vorticity at pi 

2. Integrate vorticity to find (predictor) 

3. Determine vorticity at p l+ \ (corrector) 

4. Find location p l+1 of minimum pressure in 
plane perpendicular to u; f * + i at p 1+ i 

5. Calculate quantities of interest at 

6 . i < — z T 1 

The calculation of the vortex skeleton proceeds both forward and backward from the seed 
point. The details for continuing the calculation from one point to the next are indicated 
by the numbered items in the pseudocode. Processes 1 and 2, which are shown in Fig. 2(a), 
represent the predictor stage of the algorithm. The corrector stage is summarized by pro- 
cesses 3 and 4, which are illustrated in Figs. 2(b) and 2(c). We iterate the corrector stage so 
that the pressure at position p t+i is the local minimum in the plane perpendicular to u5,+|. 
We restrict the movement in the corrector steps by limiting the angle between u3 l+ i and 
The resultant state, which is illustrated in Fig. 2(d), is equivalent to that in Fig. 2(a), except 
that the index i -f 1 replaces the index i . 

The effectiveness of the predictor-corrector scheme is illustrated in Fig. 3, in which data 
from the direct numerical simulations of Singer and Joslin [12] are analyzed. The transparent 
vortex tube is constructed with data from the full predictor-corrector method. Its core is 
indicated by the darker skeleton. The lighter skeleton follows the uncorrected integral curve 
of the vorticity. It is obtained by disabling the corrector phase of the scheme. The vorticity 
line deviates from the core, exits the vortex tube entirely, and wanders in the flow field. 

Although the general behavior of the predictor-corrector algorithm is reliable and robust, 
optimal performance of the technique requires careful attention to implementation details. 
The remainder of this subsection addresses issues that are important to the successful use 
of this method. 

2.2.1 Eliminating feeder vortices 

Because the vorticity near the edge of a vortex may be skewed with respect to the vorticity at 
the vortex center, the location of the pressure minimum in the plane perpendicular to the edge 
vorticity might not be in the vortex center. This potential mismatch of the pressure-minimum 
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and the vortex-center locations is rarely a problem in the vortex skeleton-line calculations 
because the predictor-corrector method inhibits large excursions from the vortex center. 
However, a seed point might be selected near the vortex edge. In this case, adjustment of 
the point location to the pressure minimum might still result in a point that is far from the 
vortex center. These situations result in small “feeder vortices” that spiral toward the vortex 
center. Examples of feeder vortices are illustrated in Fig. 4. We found that we can eliminate 
most of these feeder vortices by taking advantage of the asymmetry of the predictor-corrector 
method. In either forward or backward mode, the predictor-corrector scheme will converge 
to a vortex center. We can eliminate the feeder vortices by redefining the initial point (i.e., 
the point at which we start recording the elements on the skeleton line) to be the location 
after a fixed number (usually 5 to 10) of predictor-corrector steps. If the original seed point 
is located on the vortex core, the method will return to the seed on the backward integration 
path. However, if the original seed point is located on a feeder vortex, then the backward 
integration will miss the original seed point and the skeleton of the feeder vortex. Instead, 
the retraced skeleton line will closely follow the vortex center. Although the use of many 
integration points before the initial point will eliminate feeder vortices, this practice will also 
limit the minimum extent of any detected vortices. 

2.2.2 Numerical considerations 

Neither the pressure minimum nor the result of the vorticity-line integration is likely to 
be on a grid point; hence, we must choose an interpolation scheme to find the pressure 
and vorticity at arbitrary locations in the flow field. Our first experiments with trilinear 
interpolation were a disappointment. No local extrema could be found between two points 
connected by a straight line. To allow the local pressure minima to exist in the interior of a 
grid cell, a higher order interpolation method is necessary. Second-order interpolation would 
use three grid points: two from one side of the desired location and one from the other side. 
To reduce any bias from the interpolation, a four-point Lagrange interpolation is used in each 
of the three coordinate directions. When the outer points of the four-point method are not 
in the domain, two-point Lagrange (linear) interpolation is used. The complete interpolation 
scheme works quite well, although it requires more computer time than any other subroutine 
in the computer program. 

The vorticity integration can be performed with a variety of methods. First, we used a 
fourth-order Runge-Kutta approach. This produced satisfactory results; however, step-size 
optimization was difficult to automate. Instead, we developed a new technique whereby the 
point-to-point distance in the vorticity integration is always equal to the smallest dimension 
of the local grid cell. The new point location is found by advancing this distance in the 
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direction of the local vorticity vector. This procedure ensures that successive points will not 
be more than one grid cell apart, so that if the original calculation is well resolved, then the 
vorticity-line calculation will also be sufficiently resolved. The procedure also reduces the 
chance of wasting many calculations inside a single grid cell. 

2.2.3 Pressure-minimum corrector step 

The pressure-minimum correction scheme uses the method of steepest descent to find the 
local pressure minimum in the plane perpendicular to the vorticity vector. First, the vorticity 
vector at the candidate point is determined by interpolation from the surrounding grid points. 
Two perpendicular unit vectors / and g (Fig. 5) are determined in the plane normal to the 
vorticity vector. The smallest grid-cell dimension d is used as a local length scale to find 
two points in the plane in the / and g directions. The pressures at these two points and 
the current point are used to form the in-plane pressure-gradient vector. A sample point 
is chosen at a distance d from the current point in the direction opposite to the pressure- 
gradient vector. If the pressure at the new sample point is less than the pressure at the 
current point and the dot product of the vorticities at the sample and current points is 
greater than a specified value (0.90-0.95 works well), then the new sample point becomes 
the current location and the process is repeated. If the sample point has a greater pressure 
than the current point, then the distance d is halved and a new sample-point comparison 
is made. If the orginal point is still the point of minimum pressure after d has been halved 
twice, then one last check for a possible minimum is made in the direction of the positive 
gradient. (Rare circumstances with a great deal of symmetry can produce a pressure-gradient 
vector that is in the wrong direction.) Failure to find any points with lower pressure leads 
to acceptance of the current point as the next point on the vortex-core skeleton. 

2.3 The end of the line 

Vorticity lines extend until they intersect a domain boundary, but vortices typically begin 
and end inside the domain. Here, we discuss some techniques for stopping the integra- 
tion. One particularly clean termination occurs when the vortex cross section, which will 
be discussed in the next section, has zero area. This approach provides smooth surfaces for 
visualization with no abrupt vortex-tube cutoffs in the direction that is tangential to the 
vorticity. Although we have used this line-termination technique successfully, the method 
can fail to show connections between parts of the same vortex. For instance, if a low-intensity 
region exists between high-intensity regions of the same vortex, then the low-intensity region 
might not satisfy the criteria for a finite cross section. If both high-intensity regions have 
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finite cross sections, then the single vortex can be educed as two small disconnected vortices. 
Although the criteria that define the lateral edge of the vortex can be weakened to ensure 
that the low-intensity region has a finite cross section, the problem can potentially reemerge 
with a new vortex that is educed with the weaker cross-section criteria. Our resolution of 
this problem exploits the asymmetric nature of the predictor-correction method. 

Because the predictor-corrector method will follow the core of a vortex regardless of the 
criteria used to define the vortex cross section, the vortex skeleton line is continued even 
when the cross-sectional area of the vortex is equal to zero. The vortex of interest can 
either reintensify or dissipate, if the vortex reintensifies, then the continuation of the vortex 
skeleton line will provide a link between the two more intense regions of the vortex. This 
link can be visualized as a thread that connects the two disjoint regions, or the two regions 
can simply be rendered with the same color or texture. On the other hand, if the vortex 
dissipates, then the vortex skeleton line continuation will wander through the flow field and 
eventually either intercept a domain boundary or enter a new vortex. If a domain boundary 
is reached, then the elements of the vortex skeleton line that were computed after the cross- 
sectional area became zero are discarded. These same points are also discarded if the vortex 
skeleton line enters a new vortex with a nonzero cross-sectional area. To determine whether 
the new region of finite cross section is a continuation of the original vortex or an entirely new 
vortex, we march the predictor-corrector scheme backwards for the same number of steps 
taken since the previous region of nonzero cross section was exited. Some possible scenarios 
are illustrated in f ig. 6. In Fig. 6(a), the skeleton line leaves the first vortex tube at point 
Pi and continues for N steps until it encounters the second vortex tube at point I\. I he 
predictor-corrector scheme is then marched backwards N steps from P<2 to P3. The distance 
between points P\ and P3 is small relative to the distance between P \ and P 2 (a 10-percent 
criterion is used); hence, the link between P\ and P2 is most probably a low-intensity vortex, 
and we keep the thread between these vortex tubes. However, in Fig. 6(b) the vortex tube 
that ends at P 4 continues to dissipate, and the continuation of its skeleton line lacks clear 
direction and wanders through the flow field. The line intercepts another vortex tube at / 5 
after M steps. The predictor-corrector method is marched backwards A/ steps from P$ to 
P 6 . Initially, the reverse integration retraces the forward integration, but halfway between 
P s and Pq the two lines diverge rapidly and become uncorrelated. The distance from P 4 
to P& is a large fraction of the distance from P 4 to P5, so the algorithm concludes that the 
vortex tube intersected at P5 is different from the vortex tube that ends at P 4 . I he points 
on the vortex skeleton line that connect the two tubes are discarded, and the vortex skeleton 
is terminated. Finally, in Fig. 6(c), the continuation of the skeleton line of the vortex tube 
that ends at point P7 intersects the side of another vortex tube (shown as a wireframe) and 
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is immediately taken to the pressure minimum at Pg. The reverse integration for this case 
is along the axis of the new vortex tube away from the original vortex. The point P 9 is far 
from Py\ hence, the two vortex tubes are distinct from each other and the line connecting 
them is discarded. 

3 Putting meat on the bones 

The determination of an appropriate vortex cross section has been one of the more difficult 
practical aspects of this work. For isolated vortices, a simple pressure criterion to define the 
edge of a vortex works quite well, although the information content of the visualization is 
little more than that which is available from pressure isosurfaces. In regions where vortices 
interact, the pressure alone is inadequate. The low-pressure regions from two or more vortices 
can merge and distort the radius estimate of any single vortex. A similar problem arises if the 
vortex edge is defined in terms of vorticity magnitude. This particular difficulty is resolved 
if the angle between the vorticity vector on the skeleton line and the vorticity vector at any 
radial position is restricted. Any angle greater than 90° indicates that the fluid at the radial 
position is rotating in the direction opposite to that in the core. The vortex circulation T is 
defined as 

r = d; • tidA (2) 

where n is a unit vector normal to the surface S] dA is an area element on the surface 
S\ and S is a cross-sectional surface of the vortex. The vortex circulation increases as the 
cross-sectional area increases, provided that the angle between the vorticity vector at a given 
radial position and the vorticity vector on the skeleton line is less than 90°. This observation 
suggests the use of the 90° variation in vorticity vectors as a single vortex-edge criterion. 
Unfortunately, in a uniform shear flow, this criterion results in an infinite cross section 
because the vorticity at all radial locations is the same. In practice, we have found that the 
90° restriction works well in combination with a low-pressure criterion for the vortex edge. 
For the actual computation of the radial distance, the pressure and the vorticity must be 
sampled along a number of radial lines perpendicular to the vortex skeleton line. We sample 
at an increment that is equal to the minimum grid-cell dimension at the skeleton line. The 
dot product of the vorticity at the radial point with the vorticity on the skeleton line must be 
greater than zero for the 90°~variation criterion to hold. At each position, this dot product 
and the pressure are checked to ensure that they both satisfy the respective requirements. If 
either fails, the radial position along that line is the point of failure as determined by linear 
interpolation between the two most recent sample points. By sampling along radial lines 
that emanate from a point on the skeleton line, we implicitly require that the cross section 
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be starlike (i.e., each radial line will intersect the cross-section boundary only once.) Cross 
sections with more complicated geometries are truncated to a starlike form. 

Relative measures for the vortex edge are also possible. A relative measure defines the 
vortex edge as the position at which some quantity decreases below a specified fraction of its 
value on the skeleton line. Clearly, this measure will always give a finite cross section, and 
the cross section can become infinite in a uniform shear flow. We have found that relative 
measures are not as informative as the absolute measures. 

The characterization of the cross section can take many forms. The detailed aspects of 
the data are most nearly recovered if the radial location of the vortex edge is retained at 
many closely spaced azimuthal angles. If the radius of the cross section were sampled at 
10° increments in the azimuthal direction, then 36 radial distances and a reference vector 
at every position along the vortex skeleton line would need to be stored. For the minimal- 
storage approach that provides characteristic cross-sectional information, all of the radial 
distances are averaged to obtain a single scalar radius. No reference vector is needed if only 
the radius is used, so this simple method uses the least amount of storage. We have found 
that the radius is a very good cross-section descriptor for isolated vortices. When vortices 
begin to interact, the radius does not provide a good description of the asymmetric distortion 
experienced by the vortices. However, the first few coefficients of the Fourier series of the 
radial locations provides a convenient compromise between the simple radius approach and 
a full description of all radial locations. The series are easy to compute, easy to interpret, 
and allow a large range of cross-sectional shapes. In our work, we keep the constant term 
(which is the radius), the first and second sine and cosine coefficients, and a reference vector. 
Most of the cases that we have checked have a factor-of-10 drop in the magnitude of the 
first and second coefficients, a fact that suggests that the series is convergent. Figure 7 
illustrates a single cross section of a vortex educed from direct numerical simulation data. 
The shaded region represents the interior of the vortex tube, the boundary of which was 
determined by sampling along lines that radiate from the cross (which marks the skeleton- 
line location) at 1° intervals. The dashed line is a circle centered about the skeleton-line 
location with the averaged radius of the vortex tube. The solid line is the two-term Fourier 
series representation of the vortex cross section. Note how the vortex boundary represented 
by the Fourier series more nearly shows both the eccentricity and the flattening of the vortex 
tube compared to the circular cross section. 
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4 Conclusions 


The innovative use of a two-step predictor-corrector algorithm has been introduced to educe 
vortices from flow-field data. Rather than relying on any single quantity to determine the 
vortex skeleton line, the new method uses vorticity to predict the new location of the next 
element along the line and then corrects this location by using the pressure field in the 
plane perpendicular to the vorticity vector. The prediction stage resembles a portion of 
a vorticity-line calculation; the correction stage maintains the line near the vortex center. 
Unlike other approaches, our method is able to treat the skeleton line through the vortex 
core as an attractor in the flow field. 

To make the method work in all but the simplest flow fifelds, numerous side issues, 
which are common to all techniques that grow a skeleton line from seed points, must be 
addressed. This paper discusses a number of novel approaches that we have developed to 
deal with matters such as seed-point selection, feeder-vortex elimination, vortex skeleton- 
line termination, and vortex cross-section description. Sample extractions of vortices from 
various flow fields illustrate the different aspects of the technique. 

The predictor-corrector technique presented here can identify vortices in flow fields that 
are far more complex than those that we have used here for illustration purposes. An 
extension of the direct numerical simulation of Singer and Joslin [12] provides one of the more 
interesting flow configurations that we have considered. The predictor-corrector algorithm 
educes a complex tangle of many interconnected vortices that are visualized in Fig. 8. All of 
the implementation features discussed above are used to identify the vortices in this messy 
flow field. 

Although many modifications can be made to the basic scheme, the principle of using a 
vector field to predict the location of the next point and a scalar field to correct this position 
distinguishes this method from others. 
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Figure 1: Multiple realizations of same vortex tube from different seed points. 
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Figure 2: Schematic of predictor-corrector algorithm, (a) Vorticity u5 t is used to predict 
candidate point p l+1 . (b) Pressure is determined in plane perpendicular to vorticity i. 
(c) New point on skeleton line 7^+1 is corrected to pressure minimum in plane perpendicular 
to vorticity uj t +\. (d) Repeat at point p t +\ . 









Figure 3: Standard vorticity lines are compared with results from predictor-corrector scheme. 
Both results use same flow data and same computer code. Standard vorticity lines are 
obtained by disabling corrector portion of program. Transparent tube is constructed from 
data from predictor-corrector method. 



Figure 4: Feeder vortices merge with a large-scale hairpin vortex. Points that satisfy seed- 
point criteria exist on edge of vortex tube. They curve inward toward centerline and then 
follow main skeleton line. 
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Figure 5: Position correction using local pressure field. Orthogonal unit vectors / and g are 
in plane perpendicular to vorticity vector. Method of steepest descent is used to locate local 
pressure minimum. 


Figure 6(a). 




Figure 6(b). 



Figure 6(e). 


Figure 6: Schematic of line-termination procedure, (a) Forward integration from P\ to P2 
gives approximately same path as reverse integration from P 2 to P 3 , so points Pi and P2 
are connected by weak vortex, (b) Forward integration from P4 to P5 differs markedly from 
reverse integration from P$ to Pg, so vortex terminates at P 4 . (c) Forward integration from 
P7 to Pg intersects the side of new vortex tube (shown as wireframe). Reverse integration 
from Pg to Pg goes along axis of new vortex away from original tube. Distance from P 7 to 
Pg is large, so two vortices are not linked. 
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Figure 7: Comparison of different ways to represent vortex cross section. Central cross in- 
dicates point on skeleton line. Finely sampled vortex cross section is represented by shaded 
region. Dashed line is circle with averaged radius. Solid line is two-term Fourier representa- 
tion of vortex cross section. 
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Figure 8: Many interacting vortices are educed from direct numerical simulation data with 
predictor-corrector algorithm. 
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